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Abstract 

A crisp survey is given of chemical reaction networks from the perspective of general nonlinear network dynamics, 
in particular of consensus dynamics. It is shown how hy starting from the complex-balanced assumption the reaction 
dynamics governed by mass action kinetics can be rewritten into a form which allows for a very simple derivation 
of a number of key results in chemical reaction network theory, and which directly relates to the thermodynamics 
of the system. Central in this formulation is the definition of a balanced Laplacian matrix on the graph of chemical 
complexes together with a resulting fundamental inequality. This directly leads to the characterization of the set of 
equilibria and their stability. Both the form of the dynamics and the deduced dynamical behavior are very similar 
to consensus dynamics, and provide additional insights and perspectives to the latter. The assumption of complex¬ 
balancedness is revisited from the point of view of Kirchhoff’s Matrix Tree theorem. Finally, using the classical idea 
of extending the graph of chemical complexes by an extra ’zero’ complex, a complete steady-state stability analysis 
of mass action kinetics reaction networks with constant inflows and mass action outflows is given. This provides 
a unified framework for structure-preserving model reduction, and for the control (see already 1391 ) and ’reverse 
engineering’ of (bio-)chemical reaction networks. 


I. Introduction 

Network dynamics has been the subject of intensive research in recent years due to the ubiquity of large-scale 
networks in various application areas. While many advances have been made in the analysis of linear network 
dynamics, the study of nonlinear network dynamics still poses many challenges, especially in the presence of in- 
and outflows. 

In this paper, we revisit the analysis of chemical reaction networks as a prime example of nonlinear network 
dynamics, playing an important role in systems biology, (bio-)chemical engineering, and the emerging held of 
synthetic biology. Apart from being large-scale (typical reaction networks in living cells involve several hundreds 
of chemical species and reactions) a characteristic feature of chemical reaction network dynamics is their intrinsic 
nonlinearity. In fact, mass action kinetics, the most basic way to model reaction rates, leads to polynomial differential 
equations. On top of this, chemical reaction networks, in particular in a bio-chemical context, usually have inflows 
and outflows. 

The foundations of the structural theory of (isothermal) chemical reaction networks (CRNs) were laid in a series 
of seminal papers by Horn, Jackson, and Feinberg in the 1970s. The basic starting point of e.g. ED, Eg), ini 
is the identification of a graph structure for CRNs by defining the chemical complexes, i.e., the combination of 
chemical species appearing on the left-hand (substrate) and right-hand (product) sides of every reaction, as the 
vertices of a graph and the reactions as its edges. This enables the formulation of the dynamics of the reaction 
network as a dynamical system on the graph of complexes. Furthermore, in these papers the philosophy was put 
forward of delineating, by means of structural conditions on the graph, a large class of reaction networks exhibiting 
the same type of dynamics, irrespective of the precise values of the (often unknown or uncertain) reaction constants. 
This ’normal’ dynamics is characterized by the property that for every initial condition of the concentrations there 
exists a unique positive equilibrium to which the system will converge. Other dynamics, such as multi-stability or 
presence of oscillations, can therefore only occur within reaction networks that are violating these conditions. The 
main sufficient structural conditions are known as the Deficiency Zero and Deficiency One theorems, see e.g. Il20l . 
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m. For an overview of results on CRNs, and current research in this direction including the global persistence 
conjecture, we refer to lH and the references quoted therein. An important step in extending the framework of 
CRNs towards/eec/fcack stabilization has been made in ll^ : also setting the stage for further regulation questions. 

The contribution of the present paper is two-fold. First, the formulation and analysis of mass action kinetics 
chemical reaction networks is revisited from the point of view of consensus dynamics and its nonlinear versions 
0, Eol, ill. The consideration of concepts from algebraic graph theory, such as the systematic use of weighted 
Laplacian matrices, provides a framework for (re-)proving many of the previously obtained results on CRNs in a 
much simpler and insightful manner. In particular, in our previous work ll43l . we have shown how under the 
assumption of existence of a detailed-balanced equilibrium, or the weaker assumption of existence of a complex- 
balanced equilibrium (a concept dating back to Horn & Jackson EU), the weights of the graph of complexes 
can be redefined in such a way that the resulting Laplacian matrix becomes symmetric (detailed-balanced case) 
or balanced (complex-balanced case). As a result, the characterization of the set of positive equilibria and their 
stability as originating in ED, EH, lEl follows in a simple way. Moreover, this formulation allows for a direct 
port-Hamiltonian interpretation ll44l . merging CRNs with the geometric thermodynamical theory of Oster & Perelson 
E3], Gl, and leads to new developments such as a theory of structure-preserving model reduction of chemical 
reaction networks, based on Kron reduction of the Laplacian matrix ll43l . 1^ . llTSl . 

Interestingly, the nonlinearity of the chemical reaction dynamics is closely related to the fact that apart from the 
graph of complexes, another construction comes into play, namely the representation of the graph of complexes 
into the space of concentration vectors. Only if this map is the identity (corresponding to networks with complexes 
consisting of single species), the chemical reaction dynamics is linear, and, under the assumption of complex¬ 
balancedness, reduces to standard linear consensus dynamics. 

As indicated above, our approach is based on the assumption of existence of a complex-balanced equilibrium, 
generalizing the classical notion of a detailed-balanced equilibrium. Based on m a necessary and sufficient 
condition is discussed for the existence of a complex-balanced equilibrium based on the Matrix Tree theorem 
(a theorem going back to the work of Kirchhoff on electrical circuits), which extends the classical Wegscheider 
conditions for existence of a detailed-balanced equilibrium. We also make a connection with the property of mass 
conservation. Furthermore, we discuss how these results can be ’dualized’ to consensus dynamics, providing new 
insights. 

The second main contribution of the present paper is the dynamical analysis of chemical reaction networks with 
inflows and outflows. The extension of the stability theory of equilibria for reaction networks without inflows and 
outflows (called closed reaction networks in the sequel) to that of steady states for reaction networks with inflows 
and outflows (called open networks) is far from easy, due to the intrinsic nonlinearity of the reaction dynamics. 
Recently, there has been a surge of interest in open CRNs; we mention 0, 0, 0, im, ED, ED- In the present 
paper we analyze open reaction networks by revisiting the classical idea of extending the graph of complexes by 
a ‘zero’ complex ED, EQI- We will show how in this way the results based on complex-balancedness for closed 
CRNs can be fully extended to CRNs with constant inflows and mass action kinetics outflow^ In particular, while 
in Il27l . Il20l the specific properties of the zero complex do not play any role in the analysis, the present paper spells 
out the equivalence of steady states with equilibria of the extended network, and shows how due to the presence of 
inflows and outflows the set of steady states may shrink to a unique steady state, while furthermore the presence 
of steady states on the boundary of the positive orthant can be precluded. Moreover, it allows to extend the model 
reduction techniques of ll43l . 1^ . 051 to CRNs with constant inflows and mass action kinetics outflows. The 
obtained stability analysis of steady states of open CRNs with constant inflows and mass action kinetics outflows 
is one of the, up to now rare, cases of a rather complete steady state analysis of nonlinear network dynamics with 
external inputs. From a control perspective the steady state analysis of open CRNs opens the possibility of applying 
the internal model principle (see e.g. HI) to achieve output regulation for such systems with constant reference 
signals using proportional-integral controllers, for example, in the control of CSTR or gene-regulatory networks as 
in in. 

^Recently also in El the idea of adding a zero complex was followed up; however in the different context of preclusion of multi-stability 
in open CRNs. 

^This class of open reaction networks is motivated by several examples of biochemical reaction network models (see for example, the 
biochemical model in El and the examples mentioned in I2Q1 \ and also derives from the assumption that some of the complexes or species 
involved in the reactions are kept at a constant concentration (see e.g. El. El)- 
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The structure of the paper is as follows. In Section 2, based on algebraic graph-theoretical tools explored in 
ms, Ga, we give a crisp overview of the theory of closed reaction network dynamics satisfying the complex- 
balanced assumption, based on rewriting the dynamics in terms of a balanced Laplacian matrix directly linked to 
its port-Hamiltonian formulation. We provide a new perspective on the characterization of complex-balancedness 
by the use of Kirchhoff’s Matrix Tree theorem, and establish a connection to mass conservation. Furthermore, we 
indicate how ideas from reaction network dynamics may be applied to the context of consensus dynamics. Section 
3 deals with the detailed modeling of reaction networks having constant inflows and mass action kinetics outflows 
by extending the graph with an extra zero complex. Section 4 shows how the assumption of complex-balancedness 
can be extended to this case, and how this allows to derive precise results on the structure and stability of steady 
states. Section 5 provides a brief introduction to structure-preserving model reduction of open chemical reaction 
networks, based on Kron reduction of the graph of complexes. Conclusions follow in Section 6, while the Appendix 
describes how the situation of detailed-balanced chemical reaction networks can be understood as a special case of 
the complex-balanced case. 

Notation: The space of n-dimensional real vectors consisting of all strictly positive entries is denoted by R" and 
the space of n-dimensional real vectors consisting of all nonnegative entries by M". The mapping Ln : R" — 
R", X I—7- Ln(a:), is defined as the mapping whose i-th component is given as (Ln(a:))j ln(a:i). Similarly, 
Exp : R" —>■ R” is the mapping whose i-th component is given as (Exp(x))j := exp(a:i). Furthermore, for two 
vectors x,y G R” we let ^ denote the vector in R” with i-th component Finally denotes the n-dimensional 
vector with all entries equal to 1, and 0„ the n-dimensional vector with all entries equal to zero, while /„ is the 
n X n identity matrix. 

Some graph-theoretic notions (see e.g. lb)): A directed grapljl Q with c vertices and r edges is characterized by a 
c X r incidence matrix, denoted by D. Each column of D corresponds to an edge of the graph, and contains exactly 
one element 1 at the position of the head vertex of this edge and exactly one —1 at the position of its tail vertex; 
all other elements are zero. Clearly, I'^D = 0. The graph is connected if any vertex can be reached from any other 
vertex by following a sequence of edges; direction not taken into account. It holds that lankD = c — I, where I is 
the number of connected components of the graph. In particular, Q is connected if and only if ker = span 1. The 
graph is strongly connected if any vertex can be reached from any other vertex, following a sequence of directed 
edges. A subgraph of (y is a directed graph whose vertex and edge set are subsets of the vertex and edge set of Q. 
A graph is acyclic (does not contain cycles) if ker D = 0. A spanning tree of a directed graph ^ is a connected, 
acyclic subgraph of Q that spans all vertices of Q. 

II. Closed chemical reaction networks as dynamics on graphs 
A. The complex graph formulation 

Consider a chemical reaction network with m chemical species (metabolites) with concentrations x G R+, among 
which r chemical reactions take place. The graph-theoretic formulation, starting with the work of Horn, Jackson 
and Feinberg in the 1970s, is to associate to each complex (substrate as well as product) of the reaction network a 
vertex of a graph, while each reaction from substrate to product complex corresponds to a directed edge (with tail 
vertex the substrate and head vertex the product complex). 

Let c be the total number of complexes involved in the reaction network, then the resulting directed graph (y with 
c vertices and r edges is called the graph of complexe]^ and is defined by its c x r incidence matrix D. Since each 
of the c complexes is a combination of the m chemical species we define the m x c matrix Z with non-negative 
integer elements expressing the composition of the complexes in terms of the chemical species. The fc-th column 
of Z denotes the composition of the fc-th complex, and the matrix Z is called the complex composition matri^ 
It can be immediately verified that ZD equals the standard stoichiometric matrix S. The mapping Z ^ R™ 
defines a representation ll^ of the graph of complexes Q into the space R™ of chemical species (the a-th vertex 
is mapped to the a-th column of Z in R'”). Compared with other network dynamics the presence of the matrix Z 

^Sometimes called a multigraph since we allow for multiple edges between vertices. 

^In the literature sometimes also referred to as reaction graphs. 

^In 1431 . 1361 the matrix Z was called the ’complex stoichiometric matrix’. 
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constitutes a major, and non-trivial, difference; especially in case Z is not injective. The special case of Z being 
the identity matrix corresponds to single-species substrate and product reaction networks (SS reaction networks). 
The dynamics of the reaction network takes the form 

X = Sv{x) = ZDv{x) (1) 

where v{x) is the vector of reaction rates. The most basic way to dehne v{x) is mass action kinetics. For example, 
for the reaction Xi + 2 X 2 X^ the mass action kinetics reaction rate is given as v(x) = kxix\ with fc > 0 a 
reaction constant. In general, for a single reaction with substrate complex S specihed by its corresponding column 
Zs = [Zsi ■ ■ ■ Zsm\ of the complex composition matrix Z, the mass action kinetics reaction rate is given by 

. .yZsm 

ftU/i 0/2 ’ 

which can be rewritten as fcexp(ZjLna;). Hence the reaction rates of the total reaction network are given by 

Vj (x) = kj exp(Zj. Ln a:), j = 1, • • • , r, 

where Sj is the substrate complex of the j-th reaction with reaction constant k, > 0. This yields the following 
compact description of the total mass action kinetics rate vector v{x). Dehne the r x c matrix K as the matrix 
whose (j, cr)-th element equals kj, if the cr-th complex is the substrate complex for the j-th reaction, and zero 
otherwise. We will call K the outgoing co-incidence matrix (since the cr-th column of K specihes the weighted 
outgoing edges from vertex a). Then 

v(x) = KExp (Z^Lnx), (2) 

and the dynamics of the mass action kinetics reaction takes the form 

X = ZDKExp {Z'^Ln x) (3) 

The same expression (in less explicit form) was already obtained in 

It can be verihed that the cx c matrix L := —DK has nonnegative diagonal elements and nonpositive off-diagonal 
elements. Moreover, since = 0 also = 0, i.e., the column sums of L are all zero. Hence L dehnes (a 
transposed version of) a weighted Laplacian mafnhlS From now on we will simply call L = —DK the Laplacian 
matrix of the graph of complexes Q. 

B. Analysis of complex-balanced reaction network dynamics 

A chemical reaction network Q is called complex-balanced ll^ if there exists an equilibrium x* G M™, called 
a complex-balanced equilibrium, satisfying 

Dv{x*) = -LExp (Z'^Ln (x*)) = 0 (4) 

Chemically (IHi means that at the complex-balanced equilibrium x* not only the chemical species but also the 
complexes remain constant; i.e., for each complex the total inflow (from the other complexes) equals the total 
outflow (to the other complexes). Defining now the diagonal matrix 

S(x*) :=diag(exp(ZfLn(x*))).^^_ (5) 

the dynamics © can be rewritten into the form 

X = —ZC(x*)Exp (Z^Ln (-^)), C(x*) := LE(x*), (6) 

X* 

where, since Exp(Z^Ln(^)) = Ic, the transformec0 Laplacian matrix C(x*) satisfies 

£(x*)lc = 0, tlC{x*)=Q (7) 

^In ( 7 ) such a matrix L was called an out-degree Laplacian matrix. 

^In the special case imDnkerZ = {0} {deficiency zero in the terminology of 1171 ) complex-balancedness is equivalent to the existence of 
a positive equilibrium of E). 

^As shown in osi the matrix C{x*) is in fact independent of the choice of the complex-balanced equilibrium x* up to a multiplicative factor 
for every connected component of Q. 
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Hence C{x*) is a balanced Laplacian matrix (column and row sums are zero). 

Remark 2.1: Under the stronger detailed-balanced assumption Ea, m, na the Laplacian matrix C{x*) is not 
only balanced, but in fact is symmetric. In the Appendix it is discussed how the detailed-balanced situation can be 
understood as a special case of the complex-balanced one. 

Remark 2.2: Note that the vector Exp (Z^Ln (cc*)) corresponding to a complex-balanced equilibrium x* G H.™ 
defines a vector in that is in the kernel of the Laplacian matrix L. It thus follows ll22l Lemma 3.2.9] that 
the connected components of the graph (y of a complex-balanced reaction network are strongly connectec^ This 
follows as well from the fact that a graph with balanced Laplacian matrix is strongly connected if and only if it is 
connected ll2^ . The property also follows from Kirchhoff’s Matrix Tree theorem to be discussed later on. 

As shown in (generalizing the detailed-balanced scenario of ||43) a number of key properties of the reaction 
network dynamics can be derived in an insightful and easy way from the following fundamental fact. It is well- 
known a that balancedness of C{x*) is equivalent to C{x*) C'^{x*) being positive semi-definite, i.e., Ca > 0 

for all a € Based on convexity of the exponential function we can establish the following stronger property 

JM). 

Proposition 2.3: 7 ^£(a;*)Exp ( 7 ) > 0 for any 7 G R’’, with equality if and only if = 0. 

This result leads to a direct proof of a number of key properties of (| 6 ]l, which are known within CRN theory but 
proven by tedious derivations. The first property which directly follows from Proposition I2.31 l is that all positive 
equilibria are in fact complex-balanced equilibria, and that given one complex-balanced equilibrium x* the set of 
all positive equilibria is given by 


£ := {x** G R![* I S^Ln (a;**) = S''^Ln (a;*)} ( 8 ) 

In particular, the set of positive equilibria is a smooth manifold of dimension m — rank S. Another property of 
£ can be seen to be implied by the extra assumption of mass conservation, which is defined as follows. 

Definition 2.4: A reaction network with complex composition matrix Z and incidence matrix D is said to satisfy 
mass conservation if there exists p G R™ such that 

Z'^pGkevD'^, (9) 


or, equivalently, S'^/a = = 0 . 

Remark 2.5: In case Q be connected ker = span 1, and the definition of mass conservation reduces to the 

existence of p, G such that = 1. The vector /a specifies a vector of mass assignments (/a^ specifies the 
mass associated to the i-the chemical species), and the condition Z"^p = 1 means that all complexes have identical 
mass. For the general case this holds on any connected component of Q. 

Proposition 2.6: Consider a chemical reaction network as before. Then the origin 0 is on the boundary of £ if 
and only if the chemical reaction network satisfies mass conservation. 

Proof: There exists a x** G R+ with all entries arbitrarily close to 0 with S'’^Ln(a;**) = S'’^Ln(a;*) if and 
only if there exists a vector z with all entries arbitrarily close to —00 such that S’^z = S'^Ln(a;*). This, in turn, 
holds if and only if there exists a positive vector p G kerS'^ = kerZI^Z^, or, equivalently, Z"^p G kerD^. ■ 

It directly follows from the structure of the Laplacian matrix L, see e.g. 1^ . ll^ . that the dynamics Q leaves 
the positive orthant R™ invariant. Hence concentrations of chemical species remain positive for all future times. 
On the other hand, the possibility that the solution trajectories of ([^l will approach the boundary of the positive 
orthant for t —>■ 00 is not easily excluded. The reaction network is called persisten{^ if for every xq G the 
w-limit set uj(xo) of the dynamics (| 2 ^ does not intersect the boundary of R™. 

Using a result from 03, there exists for any initial condition a:o G R™ a unique x** G £ such that x**—xo G im S. 
By using Proposition 12.31 in conjunction with the Lyapunov function 

G{x) = x'^Ln 4- {x** - x)'^ (10) 

^Strong connectedness of the connected components is in CRN literature often referred to as weak reversibility (261 

^^It is generally believed that most reaction networks ai'e persistent. However, up to now this persistence conjecture has been only proved in 
special cases (cf. (D, (D, 0 and the references quoted in there). 
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it follows that the vector of concentrations x{t) starting from xq will converge to x** if the reaction network is 
persistent. The chemical interpretation is that G is (up to a constant) the Gibbs’ free energy ll34l . ll^ . 1431 '). with 
gradient vector ^{x) = Ln being the vector of chemical potentials. Hence (|6l) can be rewritten as 

BG 

i:=-Z£(a;*)Exp(Z^ —(x)) (11) 

ax 

and the‘driving forces’ of the reactions are seen to be determined by the complex thermodynamical affinitie¥^ 
'y(x) := Z'^^(x) = Z^Ln(^). Furthermore, by Proposition 12.31 equilibrium arises whenever the components 
of 7 (x) reach ‘consensus’ on every connected component of the graph of complexes Q. 


C. Port-Hamiltonian formulation 

The formulation (fTTI) admits a direct port-Hamiltonian interpretation (see e.g. El, El, El for an introduction 
to port-Hamiltonian systems). Indeed, consider the auxiliary port-Hamiltonian system 


X = Zf 


( 12 ) 


with inputs / G and outputs e G H.'^, and Hamiltonian given by the Gibbs’ free energy G defined in (fTOl) . It 
follows from Proposition 12.31 that 

/ =-£(x*)Exp(e) (13) 


defines a true energy-dissipating relation, that is, e^/ < 0 for all e G and / G satisfying (HI. By substituting 
(in]) into (El) one recovers the chemical reaction dynamics (fTTI ). 

It should be noted that the energy-dissipating relation (fTTI ) is intrinsically nonlinear, and generally cannot be 
integrated to a relation of the form / = — ^(e) for some (Rayleigh) function i? : —)> R, since the Poincare 

integrability conditions are not satisfied (unless Z is e.g. the identity matrix; see the SS reaction networks discussed 
later on). 


D. Characterization of complex-balancedness 

Complex-balancedness can be characterized as follows, cf. If45l for further details. By the definition of Ln : R!J) ^ 
R™ the existence of a complex-balanced equilibrium x* G R™, that is, LExp (Z^Ln (x*)) = 0, is equivalent to 
the existence of a vector p* G R™ such that 

LExp = 0, (14) 

or equivalently, Exp [Z'^p*) G ker L. Furthermore, we note that Exp {Z"^p*) G R(j_. 

First assume that the graph Q is connected. Then the kernel of L is 1-dimensional, and a vector p G R^_ with 
p G kerL can be computed by Kirchhojf’s Matrix Tree theorer^^ which can be summarized as follows. Denote 
the (i, j)-th cofactor of L by Cy = (—where Mi^j is the determinant of the (1, j)-th minor of L, which 
is the matrix obtained from L by deleting its i-th row and j-th column. Define the adjoint matrix adj(L) as the 
matrix with (i,_))-th element given by Gp. It is well-known that L ■ adj(L) = (detL)/c, and since deti = 0 
this implies L ■ adj(L) = 0. Since = 0 the sum of the rows of L is zero, and hence by the properties of 
the determinant it is easily seen that Gp does not depend on i; implying that Cp = pj, j = 1, ■ ■ ■ ,c. Hence the 
rows of adj(L) are given as the row vectors Pjt"^, j = 1, • • • ,c, and by defining p := (pi, ■ ■ ■ ,pc), it follows 
that Lp = 0. Furthermore, Kirchhoff’s Matrix Tree theorem says (cf. El, Theorem 14 on p.58) that Gp = pi is 
equal to the sum of the products of weights of all the spanning trees of Q directed towards vertex i. In particular, it 
follows that pk >0,fc=l,..- ,c. Moreover, since for every vertex i there exists at least one spanning tree directed 
towards i if and only if the graph is strongly connected, p G R+ if and only if the graph is strongly connected. 

Example 2.7: Consider the cyclic reaction network 


"See e.g. (H, ED, O for further information. 

^^This theorem goes back to the classical work of Kirchhoff on resistive electrical circuits 1281 : see O for a succinct treatment. Nice accounts 
of the Matrix Tree theorem in the context of chemical reaction networks can be found in ED, (S). 
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fc 2 


in the three (unspecified) complexes Ci, ( 72 , Ca. The Laplacian matrix is given as 


L = 


ki + ke 
-ki 
-ke 


-k2 
k2 + ke 
-ke 


-ke 
—ki 
ki + ke 


By Kirchhoff’s Matrix Tree theorem the corresponding vector p satisfying Lp = 0 is given as 


P = 


keke + ^2^5 + ^2^4 
kike + kiki + kike , 
kike + keke + k2ke 


where each term corresponds to one of the three weighted spanning trees pointed towards the three vertices. 

In case the graph Q is not connected the same analysis can be performed on any of its connected components. 

Remark 2.8: The existence (not the explicit construction) of p already follows from the Perron-Frobenius theorem 
EEi . ll^ Lemma V.2]; exploiting the fact that the off-diagonal elements of —L := DK are all nonnegativ^H 
Returning to the existence of p* € satisfying LExp(Z^p*) = 0 this implies the following. Let Qj, j = 
1 , • • ■ ,£, be the connected components of the graph of complexes G- For each connected component, define the 
vectors p^, - ■ ■ , p^ as above by Kirchhoff’s Matrix Tree theorem (i.e., as cofactors of L or as sums of products of 
weights along spanning trees). Then define the total vector p as the stacked column vector p := col(p^, • • • ,p^). 
Partition correspondingly the composition matrix Z as Z = [Zi • • • Zf]. Then there exists p* € R'" satisfying 
LExp p*) = 0 if and only if each connected component is strongly connected and on each connected component 


Exp (Zj p*) = Pjp ^, j = 1 ,.. ■ 


(15) 


for some positive constants I3j,j = 1, • • •£. This in turn is equivalent to strong connectedness of each connected 
component of G and the existence of constants /?' such that 

Zj/i* =Lnp^+/3'lL, j = (16) 


Furthermore, this is equivalent to strong connectedness of each connected component, and 


Finally, (fTTI) is equivalent to 


Ln p G im Z^ + ker 
D^Lnp G imD'^Z'^ = imS''^ 


(17) 

(18) 


Summarizing we have obtained 

Theorem 2.9: The reaction network dynamics x = —ZLExp (Z^Ln (a;)) on the graph of complexes G is 
complex-balanced if and only if each connected component of G is strongly connected (or, equivalently, p G M(j_) 
and (fTSl) is satisfied, where the coefficients of the sub-vectors fp of p are obtained by Kirchhoff’s Matrix Tree 
theorem for each j-th connected component of G- 

Remark 2.10: The easiest way to compute the elements pk, = 1, • • • , c, of p is by taking the determinant of 
the matrix obtained from L by deleting its fc-th row and fc-th column. 

Clearly, if im Z^ = R'^, or equivalently ker Z = 0 , then (fTTI) is satisfied for an}[3 p. 

Corollary 2.11: The reaction network dynamics x = —ZLExp (Z^Ln (a;)) is complex-balanced if and only if 
p G R(j_ and 

pr-pr-'-Pc” = 1 , (19) 


’^This implies that there exists a real number a such that —L + al^ is a matrix with all elements nonnegative. Since the set of eigenvectors 
of —L and —L -f o/m are the same, and moreover by = 0 there cannot exist a positive eigenvector of —L corresponding to a non-zero 
eigenvalue, the application of PeiTon-Frobenius to —L + aim yields the result; see oi Lemma V.2] for details. 

*‘^This is not surprising since ker Z = 0 implies zero-deficiency. 
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for all vectors a = co1((Ti , <72 , • • • , Uc) G ker Z n im D. 

Proof: Lnp G imZ^ + kerZ?^ if and only if = 0 for all a € (imZ'^ + kerD^)-*- = kerZ n imD, 

or equivalently 

0 = (Ti In Pi H-h CTc In Pc = In Pi ^ H-f In Pc° = ln(pi^ '' ■ Pc°) 

for all a G ker Z Cl ini D. ■ 

Remark 2.12: Note that the assumption of mass conservation (Definition 12.41 i may interfere with condition (fTTl i. 
Indeed, mass conservation implies ker C im (unless Z = 0), in which case (Ell reduces to Ln p G im 
In the Appendix we will indicate how the constructive conditions for the existence of a complex-balanced 
equilibrium as obtained in Theorem l2.9l relate to the classical Wegscheider conditions for the existence of a detailed- 
balanced equilibrium. 

E. SS reaction networks 

Reaction networks with single-species substrate and product complexes {SS reaction networks) correspond to 
c = m and Z = Im, in which case the dynamics (O reduces to the linear dynamics 

X = DKx (20) 

An SS reaction network is complex-balanced if there exists a positive equilibrium x* G such that DKx* = 0, 
and hence can be rewritten as 

x = -C{x*)^, C{x*):=-DKE{x*), E{x*) := diag{xl,-■ ■ ,x*m), (21) 

X* 

where C{x*) is a balanced Laplacian matrix. The set of positive equilibria of a complex-balanced SS reaction 
network is given as £ = {x** G | ZJ^Ln (x**) = D^Ln (x*)}, and thus, in case the graph Q is connected, as 
£ = {x** I X** = px*,p > 0}. 

In Remark 12.21 we already mentioned that the existence of a complex-balanced equilibrium implies that the 
connected components of the graph are strongly connected. For SS reaction networks also the converse holds, 
as follows from the above discussion, either based on Kirchhoff’s Matrix Tree theorem or on Perron-Frobenius 

theorenE- 

E Relation with consensus dynamics 

The dynamics x = —Lx = DKx with = 0 as occurring in SS reaction networks can be regarded as ‘dual’ 
to the standard consensus dynamics x = —LcX, where the Laplacian matrix Lc satisfies Lcl = 0. In a different 
context this has been explored in Q where x = —Lx with l^L = 0 was called advection dynamics. As also noted 
in Q this duality originates from a duality in the interpretation of the edges of the underlying directed graph Q. 
For X = —Lx with = 0 an edge from vertex i to j denotes ‘material flow" from vertex i to vertex j, while 
for X = —LcX with Lcl = 0 an edge from vertex i to j denotes ‘information" about vertex i available at vertex j. 
Thus in the first case the graph Q denotes a flow network, while in the latter case Q is a communication graph. 

It follows that the results described so far for flow networks can be ‘transposed’ to communication graphs and 
consensus dynamics. First of all, the Laplacian Lc with Let = 0 can be expressed as Lc = — where D 

is again the incidence matrix of Q while J (dually to the matrix K as before) can be called the incoming co¬ 
incidence matrix: the i-th column of J specifies the weighted edges incoming to vertex i. Furthermore, the idea 
of transforming the ‘out-degree’ Laplacian matrix L = —DK to a balanced Laplacian matrix C{x*) under the 
assumption of complex-balancedness of the graph (or equivalently, under the assumption of strong connectedness 
of its connected components) can be also applied to the consensus dynamics x = —LeX with Let = 0. Indeed, 
assume that the connected components of the graph Q are strongly connected. Then Kirchhoff’s Matrix Tree theorem 
provides a positive vector a G such that Lc = 0. In fact, aj is given as the sum of the products of the 
weights along directed spanning trees directed from vertex j. It follows that implying the 

another way is to make use of the result of da stating that a mass action chemical reaction network is complex balanced if it is 
strongly connected and has zero deficiency. Recall that the deficiency is defined as rankZ) — rank ZD. Since Z = 1 any SS network has 
zero-deficiency. 


conserved quantity Defining the diagonal matrix S := diag((Ti,' 


j) the transformed Laplacian 


matrix Cc := Sic is balanced, and hence C^ + Cc>0- Note that this immediately yields an easy stability proof 
of the set of equilibria £ = {x G | a; = dl, d > 0} for the consensus dynamics x = —L^x. Indeed, the positive 
function V{x) := satisfies 

_d. 
dt ' 

and thus serves as a Lyapunov function proving asymptotic stability of the set of consensus states £. Furthermore, 
for any initial condition xq the dynamics will converge to the consensus state d*!, where d* is given as d* = 
m '^7=1 with CTi, • • ■ , (Tm determined as above by Kirchhoff’s Matrix Tree theorem. 


-V{x) = —x^{L^Y, + Y,Lc)x = —x^ 


+ Cc)x < 0 , 


( 22 ) 


III. Reaction networks with constant inflows and mass action kinetics outflows 
In many cases of interest, including bio-chemical networks, reaction networks have inflows and outflows of 
chemical species. A mass action kinetics chemical reaction network with constant inflows and mass action kinetics 
outflows is described by the following extensior0 of Q 


X — ZDv{fl) T ^kD'inO'in ^d^outtlout (^) 5 ^ G 


(23) 


Here the matrices Din and Dout specify the structure of the inflows and outflows. Din is a matrix whose columns 
consist of exactly one element equal to +1 (at the row corresponding to the complex which has inflow) while the 
other elements are zero. Similarly, Dout is a matrix whose columns consist of exactly one element equal to —1 (at 
the row corresponding to the complex which has outflow) while the rest are zero. As in the closed network case, 
v(x) is the vector of (internal) mass action kinetics reaction rates given by v(x) = ATExp (Z^Ln (x)). Furthermore, 


is a vector of constant positive inflows, while Vout(x) G 


is a vector of mass action kinetics outflows 


described by mass action kinetics as 


Doutt^out(^) — ^outExp Ln(x)), 


(24) 


where Aout is a diagonal matrix with non-negative elements given by the mass action kinetics rate constants. 

A classical idea due to ED is that by the addition of an extra complex the reaction network (|2^ can be represented 
as a closed reaction network on the extended grapf0 In fact, an extra complex is added in such a way that the 
edges from the extra complex to the ordinary complexes model the inflows into the network, while the edges towards 
the extra complex model the outflows of the network. The complex composition matrix Zzero corresponding to the 
extra complex is defined to be the m-dimensional zero column vector, and the extra complex is therefore called the 
zero complex. Hence the zero complex serves as a combined ‘source and sink’ complex, which does not contribute 
to the overall mass. As a consequence, the extended network cannot satisfy mass conservation, cf. Definition |9l 
The resulting graph, consisting of the original graph of complexes together with the zero complex, is called the 
extended graph of complexes of the open reaction network (|2^ . and has complex composition matrix 

Ze = \Z .^zero] = \Z O] 

The incidence matrix of the extended complex graph, denoted by Dg, is given as 

B 


D,= 




with Dzero a row vector corresponding to the zero complex, while in the notation of 

B = \D Din Dout] 

Now define the c-dimensional column vector Dm as 


L\n — — D'u 




(25) 


(26) 


*®Note that {ID formalizes a situation of direct in- and outflow of some of the chemical complexes in the reaction n etwo rk. Modeling of in- 
or outflows of single chemical species which do not already appeal* as complexes in the graph need to be incorporated in {ID by the introduction 
of extra complexes. For other scenarios of open chemical reaction reaction networks such as continuous-stirred tank reactors with convective 
in- and outflows we refer to e.g. {ID 

’^Similar ideas of adding vertices to the graph are used in network flow theory; see e.g. (D- 
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Furthermore, let Lout be the c-dimensional column vector whose i-th element is equal to minus the i-th diagonal 
element of Aout- Then extend the cx c Laplacian matrix L of the graph of (ordinary) complexes to an (c+1) x (c+1) 
Laplacian matrix of the extended graph of complexes as 


L, := 


L + Aout 
-Tout 


Ln 

(5iTl 


(27) 


where (5in > 0 equals minus the sum of the elements of L-m- By construction Lg has non-negative diagonal elements, 
non-positive off-diagonal elements, while its columns sums are all zero. 

It follows that the dynamics (|2^ of the mass action reaction network with constant inflows and mass action 
kinetics outflows is equal to the mass action kinetics dynamics of the extended graph of complexes with extended 
stoichiometric matrix Se = ZgDe- Indeed, since ^zero = 0 


X = Z[Dv{x) + DinVin + -Doutt^out (a^)] = ZBVe{x) = ZeDeVe{x) = SeVe{x), (28) 


where 


Ve{x) 


v{x) 

Vin 

t^out (a^) 


with v{x) = FfExp (Z^Ln (a;)). Furthermore, by using (l26l l. (l24l i. dZTl i. 


X 


Z[Dv(^x') -f -Din'Uin T -Tloutfout (a^)] 


-[z 


0] Le 


Exp (Z'^Ln (x)) 
1 


—ZeLgExp (ZgLn (x)) 


(29) 


(30) 


Example 3.1: As a simple example consider a reaction network with x € consisting of one reversible 
reaction with forward and reverse reaction constants /c_|_, fc_ > 0, where there is a constant inflow fcjn towards the 
first complex Xi and mass action kinetics outflow fcouta; 2 a;| out of complex X 2 + 2 X 3 , that is 





The complex composition matrix Z for this case is given by 


Z = 


1 

0 

0 


0 

1 

2 


The Laplacian matrix of the internal reversible reaction (split into a forward and reverse reaction) is 


L = 


fc+ —k- 

—/c_i_ k_ 


Together with the zero complex this corresponds to the Laplacian of the extended graph of complexes 


Lg = 


k+ —k- —k'u 
-k+ fc_ -I- Tout 0 

0 Tout Tin 


and the following dynamics of the reaction network as in (l23t 


X = Z 


—T_|_ T_ 
T+ — T_ 


Xl 

X 2 x 1 


0 

Touta;2a;i 
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IV. Analysis of reaction networks with constant inflows and mass action kinetics outflows 

Based on the formulation of the previous section we can extend the results concerning the stability of closed 
complex-balanced reaction networks as described before to the case of reaction networks with constant inflows and 
mass action kinetics outflows. As before we note that the representation (l30l l implies that the positive orthant K.™ 
is invariant for ( |2^ . 

Definition 4.1: An x* € R™ is called a steady-state of the reaction network with constant inflows and mass 
action kinetics outflows given by (1^51) if 

Z[Dv{x*) + An'yin + -Doutt^out(a;*)] = 0 (31) 

An X* € R™ is called a complex-balanced steady-state if 

Dv{x*) Anfin + -Doutfout(T:*) = 0 (32) 

If there exists a complex-balanced steady-state x* G R™ then the open reaction network (I23|) is called complex- 
balanced . 

Note that, like in the case of closed networks, at a complex balanced steady state the total inflow from every 
complex is equal to the total outflow from it. 

The definition of a complex-balanced steady state x* can be succinctly written as Bve{x*) = 0, with B given by 
(I 25 I 1 and Ve given by (|29] |. We have the following simple but crucial observation showing that complex-balanced 
steady states for (l23T l are actually complex-balanced equilibria of the extended network, and conversely. 
Proposition 4.2: x* is a complex-balanced steady-state, i.e., Bve{x*) = 0, if and only if DeVe{x*) = 0. 

Proof: Since = 0 the last row of is dependent on its first c rows, that is, the rows of B. Hence 

Ve{x*) € keiL>e if and only if Ve{x*) € kei B. ■ 

Remark 4.3: Note that this proposition does not remain true if we would consider instead of a single zero complex 
e.g. a source and a sink complex. 

If the network with constant inflows and mass action kinetics outflows has complex-balanced steady state x* then, 
similarly to Q for closed complex-balanced reaction networks, we define the diagonal matrix 

Ee{x*) := diag(exp(ZfLn(A))).^^_ 


and rewrite 


where 


'diag (exp(Zf Ln (a:*))) O' 


■S(x*) O' 

0 ’ ’ 1 


0 1 


DeVe{x) = -jCe{x*)Exp 


Z^Ln(- 

0 


£e(x*) := L,Ee{x*) 


(i -f Aout)'^(T^*) Lin 
Lont^ix*) (5in 


(33) 


Note that Exp 


0 


= Ic+i. Hence, the existence of a complex-balanced steady state x* implies by 


Proposition 14.21 that Ce{x*)'kc+i = 0. Hence, similarly to the previous section, Ce{x*) satisfies = 0, 

as well as £e(x*)lLc+i = 0, and thus defines a balanced weighted Laplacian matrix for the extended graph of 
complexes. The fact that the sum of the elements of the last row of Ce{x*) is zero amounts to the equality 
LoutS(a;*)lLc + i^in = 0, or equivalently. 


LoutExp(^^Ln(x*)) = if An, 


(34) 


which can be interpreted as a mass-balance condition: at steady state the total inflow in the reaction network is 
equal to the total outflow. 

We obtain the following refined version of Proposition 12.31 
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Theorem 4.4: Define Ce{x*) as above. Then 


7j£e(a;*)Exp (7e) > 0 (35) 

for all 7e, while equality holds if and only if = 0. Furthermore, if 7e has last component zero, i.e., is of the 
form 


7e = 


7 

0 ’ 


(36) 


then equality holds if and only if i?^7 


= 0, or equivalently 

= 0, A^7 = 0, = 0 


(37) 


Proof: Only the last statement remains to be proved. This follows by noting that if 7e is given as in (l36l) then 
£)^7e = 0 if and only if 5^7 = 0. ■ 

We obtain the following basic theorem extending and refining the results of the previous section from closed 
networks to open reaction networks. 

Theorem 4.5: Consider a mass action kinetics reaction network with constant inflows and mass action kinetics 
outflows (l23l) . for which there exists an x* G R™ satisfying (l32l) . Then 

(1); The set of positive steady states is given as 

{x** G R™ I {x**) = B^Z'^Ln (x*)}. (38) 


and all positive steady states are complex-balanced. 

(2) : If every component of the graph of complexes is connected to the zero complex (or equivalently, if the 
extended graph of complexes is connected) then the set of steady states is given as 

{x** G R™ I {x**) = Z'^Ln (x*)}. 

In particular, if additionally Z is surjective then the steady state x* is unique. 

(3) : For every xo G M™, there exists a unique xi G f with xi — xo G Iras'. The steady state xi is locally 
asymptotically stable with respect to initial conditions xq with xi — xq G im S. Furthermore, if the network is 
persistent then xi is globally asymptotically stable with respect to all these initial conditions. 

Proof: (1): (l38l l follows from the characterization of the set of equilibria of a closed network, cf. since 
the transpose of the stoichiometric matrix for the extended network is given as Sj = Zj = Z^. That every 

positive steady-state is complex balanced can be proved similar to the case of the closed networks case. 

(2); Let x** be a positive steady state, and define 7(x**) = Z^Ln (^^). By the first part of the theorem 
this means that i3^Z^Ln(x**) = i?^Z^Ln(x*), which is the same as lP"^{x**) = 0, or equivalently dTTl i. In 
particular, Dinjix**) = 0 and Dout'lix**) = 0 , and thus the components of 7(x**) corresponding to the complexes 
directly linked to the zero complex are zero. Furthermore, since B'^j{x**) = 0, it follows that the components of 
7(x**) corresponding to each of the connected components of the extended graph are equal. Hence if the extended 
graph of complexes is connected, then 7(x**) = 0, which is the same as Z^Ln(x**) = Z^Ln(x*). In particular, 
if Z is surjective then this implies that the steady state x* is unique. 

(5): This follows directly from the closed network case. ■ 

It can be concluded from Theorem 14.51 that the presence of inflows and outflows has the tendency to ‘shrink’ 
the set of positive equilibria for the closed network to a smaller set of positive steady states; in fact, to a singleton 
if the extended graph is connected and Z is surjective. Furthermore, as shown in the following proposition, if the 
extended graph is connected, no steady states can occur at the boundary of the positive orthant R!p, implying that 
the reaction network is automatically persistent. 

Proposition 4.6: Consider a reaction network with constant inflows Uin G R^. and mass action outflows (l23l) . 
which is complex-balanced. If the extended graph of complexes is connected, then there are no steady states at the 
boundary of R™. 

Proof: Assume by contradiction that there exists a steady state Xb G R™ with at least one component (say the 
f-th one) equal to zero. Then consider a complex C containing this f-th species. Because Xbi = 0 the outflows from 
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complex C are zero, and by complex-balancedness this means that also all inflows to it are zero. From Remark |2^ 
it follows that the extended graph of complexes is strongly connected. Hence there exists a directed path of reactions 
n starting from the zero complex and ending at C. Now consider the complex which is preceding the complex 
C in this path. Then its outflows are zero, and therefore by complex-balancedness also its inflows. Repeating this 
argument this shows that along H the inflow from the zero complex is zero, which yields a contradiction. ■ 
Example 4.7: Consider the reaction network in Example 13. II with the Laplacian matrix of the extended graph of 
complexes given as 


Le 


k+ -/c_ — fcin 

-k+ fc_ + kout 0 

0 ^out ^in 


A complex-balanced steady-state x* = satisfies the equations 


k+ 

-k+ 


-k_ 

k— -f /Cout 



X *1 


/Cin 


X* 2 {xlf 


-1 

o 

_1 


or more explicitly 


koutxlixlf = kin, k+xl = kin + k_x* 2 {xl)'^ 


Hence the network is complex-balanced if ^ 0 and k^nt ^ 0 (and also in the degenerate case kin = kout = 0). 
The mass-balance condition (IMl) in this case amounts to 

kin = koutX*2{,xlf. (39) 


If kin 7^ 0 and fcout ^ 0 then the set of steady states of the network is 1-dimensional. Note on the other hand that 
the set of equilibria for the case without inflows and outflows {hm = kout = 0) is 2-dimensional; in line with the 
observation that the addition of inflows and outflows has the tendency to shrink the set of steady states as compared 
to the set of equilibria. Finally, note that the matrix ICe{x*) in this example equals diag (if , a;2(a;|)^, 1), while the 
resulting matrix Ce{x*) is given by 


Ce{.X*) 


k+xl -k_xl{xlY -kin 

-k+xl {k-+ kout)xl{xlY 0 
0 -koutxl(xlY kin 


V. Structure-preserving model reduction of open reaction networks 

As detailed in the previous section, the dynamics of a complex-balanced chemical reaction network with constant 
inflows and outflows governed by mass action kinetics (’open reaction network’) can be written in terms of the 
balanced Laplacian matrix given by ( l33l l as follows 

X = -Ze£e(a;*)Exp(ZjLn(-^)). (40) 

X* 

This specific form allows for application of the model reduction method discussed in llTSll . ll^ ; see also ||4^ for 
the detailed-balanced case. This method is inspired by the Kron reduction method of resistive electrical networks 
described in 1^ ; see also HI]. The speciality of this method is that it is structure-preserving in the sense that 
the reduced model corresponds to a complex-balanced chemical reaction network governed by mass action kinetics 
just like the original model. To make the paper self-contained, we briefly describe the method below. For a detailed 
description of the method, the reader is referred to m, 136). 

Let V denote the set of vertices of the graph of complexes. We perform model reduction by deleting certain 
complexes in the graph of complexes, resulting in a reduced graph of complexes. We ensure that the set of complexes 
that are deleted does not include the zero complex, because otherwise the reduced network corresponding to an 
open network would become a closed network. Deletion of a complex is equivalent to imposing the complex¬ 
balancing condition on it, i.e., the condition that the net inflow into the complex is equal to the net outflow from 
it. Consider a subset 120 c V of dimension c -f 1 — c that we wish to delete in order to reduce the model. Without 
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loss of generality, assume that the first c rows and columns of Ce{x*) and the first c columns of Zg correspond to 
Vr '■= V\Vo- Consider the resulting partition of Ce{x*) given by 


'C^^{x*) Ci2{x*) 

C2l{x*) C22{x*) 


(41) 


where Cii{x*) € £i 2 (x*) G ^^><(=+1-2), £21(3 :*) G K(^+i-2)xS and £22(2 :*) G IR(‘:+i-s)x(c+i-c)^ 

resulting partition of Zg given by Zg = [Zi Z 2 ]. Then write out the dynamics (l40l) as 


X = - [Zi 


^2] 


£ii{x*) 

£ 2 i{x*) 


£i 2 {x*) 

£22(3;*) 


[Exp f^fLn (^J)I 
[Exp (^2^Ln (^))J 


Let Ce{x*) denote the Schur complement of £e{x*) with respect to the indices corresponding to Vo- Consider now 
the auxiliary dynamical system 


m 


£ii(x*) 

£i2(x*y 


Wi 



C 2 l(x*) 

£22(x*) 


W2 


Note that the complex-balancing condition on the complexes in Vo can be imposed by setting the constraint j/2 = 0. 
This results in the equation 

W2 = -C22{X*)~^ C2l{x*)wi, 


leading to the reduced auxiliary dynamics defined by the Schur complement 

yi = -(£11(3:*) - Ci 2 ix*)C 22 {x*)~^ C 2 l{x*))wi = -Ce{x*)wi (42) 

Substituting wi = Exp (Z^Ln (x)) in the above equation and making use of i = Ziyi + Z 2 y 2 = we then 

obtain the reduced model given by 

X = -Zg£e(x*)Exp (Zg^Ln . (43) 

where Zg := Zi. The following proposition ensures that £e(x*) obeys all the properties of the weighted Laplacian 
matrix of a complex-balanced reaction network corresponding to a graph of complexes with vertex set Vr- 

Proposition 5.1: Consider an open complex-balanced network with dynamics given by equation (l40l i. With V, 
Vo and £e as defined above, the following properties hold; 

1) All diagonal elements of £g(x*) are positive and off-diagonal elements are nonnegative. 

2) lJCe{x*) = 0 and £g(x*)lc = 0, where c := c-f 1 — dim(Vo). 

If £ and £ denote the set of steady-states of the original and the reduced networks described by (l40l) and (03) 
respectively, then £ <Z £. 

Proof: See proofs of ll^ Propositions 5.1,5.2] ■ 

From Proposition 15. II it follows that the reduced network (03 corresponding to a complex balanced network (03 
is also complex balanced. Complexes belonging to a certain connected component remain in the same connected 
component if not deleted. Thus mass conservation is preserved under our model reduction procedure. 

Finally, we remark that for the application of our model reduction method, one can also start directly from 
the form of equations (133 given by i = ZgLgExp (ZjLn (x)), instead of the form (03 that uses the balanced 
Laplacian. Let £g denote the Schur complement of Lg with respect to the indices corresponding to Vo- Consider 
the reduced model given by 

X = ZgLgExp (ZjLn(x)) 

and note that it is the same as the reduced model (143 . 
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VI. Conclusions and outlook 


We have discussed mass action kinetics chemical reaction networks as a challenging example of large-scale and 
nonlinear network dynamics, and have pointed out similarities with (nonlinear versions of) consensus dynamics. A 
fundamental difference resides in the complex composition matrix Z, which defines a representation of the graph of 
complexes (into the space of chemical species). Kirchhoff’s Matrix Tree theorem has been discussed as an insightful 
way to compute the kernel of the Laplacian matrix, which, among others, yields an explicit characterization of the 
existence of a complex-balanced equilibrium. Also the relation to mass conservation has been pointed out. 

For a particular class of open reaction networks, namely those with constant inflows and mass action outflows, 
a detailed stability analysis has been obtained by exploiting the notion of zero complex. By using the graph- 
theoretical techniques that we have used earlier to analyze closed complex-balanced reaction networks this 
leads to a complete steady state stability analysis for this class of open reaction networks. Our results imply the 
intuitively obvious fact that the presence of inflows and outflows has the tendency to shrink the set of positive 
equilibria to a smaller set of positive steady states, and leads to the vanishing of possible steady states at the 
boundary of the positive orthant. This can be related to the feedback stabilization problem studied in 1^ . as well 
as internal model control. 

An important extension of our results concerns the consideration of other types of kinetics, in particular Michaelis- 
Menten kinetics; see already for the closed network case. Furthermore, the framework described in this paper 

can serve as a starting point for the inclusion of regulatory networks thus leading to direct control and ’reverse 
engineering’ questions. 


VII. Appendix: Detailed-balanced reaction networks 

The assumption of existence of a complex-balanced equilibrium can be strengthened to the existence of a detailed- 
balanced equilibrium. In this case we start with a directed graph of complexes T-L with c complexes and p edges, 
where each edge corresponds to a reversible reaction, cf. ms. Assuming again mass action kinetics, the reaction 
rate (x) of each j-th reversible reaction is given as the difference 

Vj(x) = exp{Zg^Lnx) — k~ exp{Z^^Lnx), 

where Sj is the substrate and Vj the product complex, and where k'^ and k~ are respectively the forward and 
reverse reaction constants of the reversible reaction. Note that v^^ix) may take positive and negative values, in 
contrast with the previously considered case of irreversible reaction rates Vj{x) > 0. A reversible reaction network 
can be brought into the irreversible form as discussed before by defining the directed graph Q as having the same 
vertex set as H but with twice as many edges: every edge {i,j) of H is split into two edges (of opposite orientation) 
(i,j) and {j,i) of g. 

A reversible mass action kinetics reaction network with graph of complexes H is called detailed-balanced if 
there exists an x* € K™ satisfying v'^{x*) = 0, i.e. 

exp [zl.hn (x*)^ = kj exp (^Z^.Ln (x*)^ , j = 1, • • • , r 

It is immediate that detailed-balancedness is a special case of complex-balancedness, with the reaction rates in the 
two opposite edges of Q corresponding to a single edge of H being equa0. (Instead of having the total sum of 
inflows to be equal to the total sum of outflows for every complex.) 

k'f _ 

Defining the equilibrium constants = -^ of each reversible reaction, and the vector := • • • , 

J Kj 

it can be shown ifTSl . BtI . Il4^ that detailed-balancedness is equivalent to the Wegscheider conditions 

LnA'®'^ € imD^Z^ = im^.^, 

with D-h the incidence matrix of the graph T-L, and S-h the stoichiometric matrix corresponding to Ti (i.e., every 
column of Su corresponds to a reversible reaction). The assumption of detailed-balancedness implies that all 
equilibria are actually detailed-balanced, and that we may define the conductances of the j-th reversible reaction as 

Kj{x*) := fc+ exp {z^.hTL {x*)^ = kJ exp (a:*)) , j = 1, •. • , r (44) 

^^Thermodynamically the assumption of detailed-balancedness is well-justified; it corresponds to microscopic reversibility ED. 
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(See ini, m for a discussion regarding the similarities of these constants with conductances in other physical 
networks.) It is readily seen that the detailed-balanced assumption is equivalent to the transformed Laplacian matrix 
C{x*) := L'E.{x*) for the graph Q being symmetric, with the (i,_))-th = (j,i)-th element being equal to the 
conductance of the reversible reaction between the i-th and the j-th complex. This means that C{x*) can be written 
as 

Cix*) = DnlC{x*)Dl^ 

wher^H lC'{x*) is the diagonal matrix of conductances Kj{x*),j = 1, • ■ • ,r. Hence ll43ll the dynamics takes the 
form 

i = -ZDnlC^x*)Dl,E^piZ^hn{^)). 

For Z = I this amounts to symmetric consensus dynamics on the graph T-L without its orientation. 

Finally we will indicate the connections of the Wegscheider conditions mentioned above to the characterization 
of complex-balancedness as obtained by the Kirchhoff’s Matrix Tree theorem. For further details we refer to ||45]| . 
Consider a complex-balanced reaction network, with Laplacian matrix L = —DK. Compute based on Kirchhoff’s 
Matrix Tree theorem the vector p € satisfying Lp = 0, leading to the transformed balanced Laplacian matrix £ 
given as £ = Ldiag (pi, • • • , pc)- In case the reaction network is detailed-balanced it follows that this transformed 
Laplacian matrix £ is actually symmetric instead of just balanced. Symmetry of £ can be seen to be equivalent 
to the weakened Wegscheider conditiono (only depending on the structure of the graph of complexes, and not on 
the composition of the complexes) 

Ln/C®'^ G im£)^ 

Example 7.1: Consider the network described in Ex. 12.71 The transformed Laplacian matrix is computed as 



ki + kg 

-k2 

-ks 


^ 3^5 + ^ 2^5 + ^2^4 

0 

0 

£ = 

-ki 

^2 + ^3 

—ki 


0 

kik^ kiki -\- kik^ 

0 


-ke 

-ks 

ki -1- /cs 


0 

0 

kiks k^ke -f k2kQ 


(fci-I-fc3)(fc3A:5-f/C2fc5-f ^2^4) — fc2(fcifc5-I-fciA:4-I-/C4/C6) —k^ikik^-\-k^kQk^k^) 

-ki{k3k^k2k^k^ki) {k2k3){kik^kik4^-\-k^k^) —k^^^kik^-\-k^k^k2k^) 

-kQ{kzk^k2k^k2k4) —k'i{kxk^kik^k^ko) {k^k^){kik'ik^k^,k2kQ) 


This matrix is symmetric if and only if 

kiksk^ = k2kikQ 

On the other hand, Ln G amounts to 


In K 


■-1 

0 

1 ■ 

ln|^ 

K 4 

G im 

1 

-1 

0 

ln|^ 

kQ_ 


0 

1 

-1 


which reduces to In -f In -f In || = 0, and hence to the same condition kik^k^ = k 2 k 4 kQ. 


*^It can be shown ED that the matrix is independent of the choice of the thermodynamic equilibrium x* up to multiplicative factor 

for every connected component of "H. 

^^As shown in ED the weakened Wegscheider conditions ai‘e also equivalent to the notion of formal balancing introduced in (H as a 
formalization of the ’circuit conditions’ of mi 
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